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Abstract 

This series of papers is devoted to the formulation and the approximation of coupling 
problems for nonlinear hyperbolic equations. The coupling across an interface in the physical 
space is formulated in term of an augmented system of partial differential equations. In 
an earlier work, this strategy allowed us to develop a regularization method based on a 
thick interface model in one space variable. In the present paper, we significantly extend 
this framework and, in addition, encompass equations in several space variables. This new 
formulation includes the coupling of several distinct conservation laws and allows for a possible 
covering in space. Our main contributions are, on one hand, the design and analysis of a 
well-balanced finite volume method on general triangulations and, on the other hand, a proof 
of convergence of this method toward entropy solutions, extending Coquel, Cockburn, and 
LeFloch's theory (restricted to a single conservation law without coupling). The core of our 
analysis is, first, the derivation of entropy inequalities as well as a discrete entropy dissipation 
estimate and, second, a proof of convergence toward the entropy solution of the coupling 
problem. 

1 Introduction 

Objective of this paper 

This is a continuation of a series of papers [111 1121 113) devoted to coupling techniques for nonlin- 
ear hyperbolic equations. In the present paper, we deal with the coupling of multi-dimensional 
hyperbolic equations, based on an arbitrary partition of the physical domain. The main moti- 
vation stems from the study of complex systems resulting from the combination of elementary 
components modeled by different equations. Indeed, each component may be subject to physical 
phenomena involving fairly different time and space scales. Tackling this multiscale problem with 
sufficient accuracy and efficiency requires to consider distinct physical models for the description 
of each component, so as to end up with a suitable description of the whole physical system. For 
instance, large-scale power plants provide a typical example of interest |29| . Describing the evolu- 
tion in time requires the exchange of transient informations at each physical boundary separating 
two distinct hyperbolic models. These transient informations or boundary conditions are referred 
hereafter to as coupling conditions. 
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This problem seems to be rather new in the applied mathematical community. Its analysis was 
initiated by Godlewski and Raviart |28j for scalar equations in one space variable. Therein, the 
coupling problem is formulated in terms of two initial boundary value problems (IBVP) supple- 
mented with coupling boundary conditions at a given (infinitely thin) interface. These boundary 
conditions are stated in such a way that in "most cases" they ensure the continuity of the main 
unknown, at least, roughly speaking, as long as no wave from the left- and right-hand prob- 
lems interact at the interface. If this condition does not hold, one says that the interface is 
resonant. In Ambroso et al. El El [5], quite general continuity conditions based on a nonlinear 
transformation of the unknown were investigated. Following earlier investigations by LeFloch and 
collaborators [23J EH EH ESI E3 EI] on undercompressive shocks and interfaces, nonconservative 
hyperbolic systems, and boundary value problems, we stress that additional information coming 
from physical modeling is necessary in order to single out the relevant continuity conditions (or 
transmission condition) at the interfaces. Various conditions were introduced and studied in a 
variety of physical frameworks, ranging from gas dynamics [3 to multiphase flows |UI1]. 

Thin interface versus thick interface 

We briefly mention some transmission conditions of interest when the coupling invoves two Euler 
systems with distinct pressure laws. Typically, one imposes the continuity of the density p, velocity 
component u, and pressure p, or else the continuity of the convervative variables (p, pu, pE) (where 
E denotes the total energy). These conditions determine the class of constant solutions in the 
time and the space variables, and have either constant density, velocity, and pressure, or else 
constant density, momentum, and total energy. In both cases, the proposed coupling conditions 
are nonconservative, since the total mass of density, momentum, and total energy do vary with 
time. A fully conservative coupling may turn relevant in some applications, as was addressed in 
[5] (following (2H1) via suitable a relaxation method. 

The resonance phenomena, likely to take place around thin interfaces, brings a main difficulty 
in the mathematical analysis of coupled initial boundary value problems. Solutions can be shown 
to exist under general conditions but resonance generally comes at the expense of uniqueness. 
We refer the reader to [9] for a discussion of scalar equation and to [3] for a distinct behavior 
exhibited for characteristic but non-resonant interfaces. A selection criterion for discontinuous 
solutions, therefore, is required. Recall that, for the fully conservative coupling, several distinct 
entropy criteria have been proposed, each selecting a distinct weak solution in agreement with the 
physical context. (See [14] for a review and [4T | l30 l 17]). 

To deal with general transmission conditions, a macroscopic selection principle analogous to the 
entropy inequalities is not available and one needs a detailed description of microscopic mechanisms 
coming with suitable regularizing procedures. In |10LI11L[T2] . we introduced an alternative modeling 
for the coupling problem associated with two hyperbolic equations in one space variable. This 
alternative method relies on the introduction of an augmented PDE (partial differential equations) 
formulation that avoids the detailed description of the interfaces. The proposed formalism is based 
on an additional unknown, the color function which takes values in the range [0, 1]. Extreme values 
and 1 are devoted to restore the left- and right- problems to be coupled, while intermediate 
values may serve to model a smooth transition from one problem to the other. 

Outline of this paper 

The interest in this augmented formulation comes from its very capability to support various 
regularization mechanisms. Viscous perturbations were introduced by the authors [111 112] for 
scalar problems and, specifically, a self-similar approach was developed, which allows for the 
study of the existence and uniqueness of solutions to the coupled Riemann problem in the limit of 
vanishing viscosity. The analysis has been carried out for a general class of systems [11] and led to 
an existence theory under fairly general assumptions. In |12j . the analysis of the internal structure 
of resonant interfaces was performed and led us to a characterization of the set of admissible 
Riemann solutions. Despite of the viscous mechanisms a failure of uniqueness may be observed 
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for resonant infinitely thin interfaces. 

Riemann solutions may be indeed understood as describing the long time asymptotic of the 
solutions of the Cauchy problem. Failure of uniqueness for thin interfaces just reflects the property 
that distinct regularizations of thin interfaces may give rise to different solutions and thus with 
a distinct long time behavior. This observation has motivated a second regularization procedure 
based on thick interfaces. 

Thick interfaces within the augmented PDE framework are based on a regularization of the 
discontinuous color function, considered in the thin regime. This approach has been introduced by 
the authors first within the framework of two coupled conservation laws in one space variable |13j . 
Existence and uniqueness for the coupled Cauchy problem was proven for general initial data with 
bounded sup-norm. One of the main ingredients of proof was the design of a well-balanced finite 
volume method. The well-balanced property means that the exact constant solutions selected by 
a given transmission condition are exactly preserved at the discrete level, whatever choice is made 
for the regularized color function. This consistency property is of central importance. 

In the present paper, we introduce a framework which covers coupling problems in several 
space variables and with distinct hyperbolic equations, allowing for possible covering in space. An 
outline of this paper is as follows. 

• In Section 2, we show how to extend the two existing coupling frameworks in one space 
variable to the coupling of two distinct hyperbolic equations in several space variables. We 
then show how to extend the augmented PDE formalism to encompass the case of several 
hyperbolic equations with possible covering. In our approach, a vector-valued color map is 
introduced so that each component is associated with one of the equations and takes values 
in the interval [0, 1]. The specific definition of the regularized color function provides us with 
a transition from an equation to another (possibly more than one). 

• We check the existence and uniqueness of entropy solutions to the coupled Cauchy problem 
(with initial data in L°°) under fairly general assumptions on the transmission conditions 
and the equations under consideration. 

• Next, in Section 3, we design a robust and flexible finite volume framework based on general 
triangulations. Importantly, by construction, the proposed method is well-balanced and our 
strategy for achieving the well-balanced property is an extension of the subcell reconstruction 
approach (analyzed by Bouchut in a different context [5]). In particular, we introduce two 
distinct meshes: the first one, the primal mesh, describes the main coupled unknown. The 
second mesh, referred to as the dual one, is built from the primal mesh and carries the 
approximation of the color function. A comprehensive derivation of this dual mesh is also 
proposed. 

• In Sections 4 and 5, we then derive a sup-norm estimate, and observe that a uniform estimate 
on the total variation seems to be out of reach, due to the subcell reconstruction procedure. 
Consequently, we propose to use DiPerna's framework based on entropy measure-valued 
solutions and, by deriving suitable entropy inequalities and entropy dissipation bounds, we 
establish the strong convergence of the proposed method. 

• Finally, in Section 6, numerical experiments are presented which concern problems with 
covering in space and, therefore, highlight the interest of the new coupling strategy. 

2 A framework for multi— dimensional coupling 
2.1 Coupling of two systems 

Pasting together two initial boundary value problems 

In this section, we introduce the coupling problem associated with two hyperbolic equations cou- 
pled at a given interface. At this stage, it suffices to think of an hyperplane, say {x\ — 0}. We will 
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extend two distinct coupling strategies that have been developed in a single space variable. The 
first procedure consists in modeling the coupling problem as two initial boundary value problems 
(IB VP) with time dependent boundary conditions prescribing the evolution of traces of the coupled 
solutions on both sides of the hyperplane {xi = 0}. In contrast, the second strategy introduced in 
[TT1 ITS"] is based on augmented PDE systems, and handles the coupling problem as an initial 
data problem written over the entire space R d . This new framework brings mathematical and 
numerical advantages, pointed out at the end of this section. 

Consider an hyperplane of M d with unit normal vector v 6 lR d , we denote W = {i£ ]R d jx.v = 
0}, partitioning Ft d into two half-domains X>_ = {x £ M d jx.v < 0} and V + = {x £ Wl d jx.v > 0}. 
In each open subdomain, a distinct conservation law is prescribed: 



d t w 



i=i 



d. r 



: (w) = 0, w(t,x)eM, t>0, x£V ± , 



(2.1) 



where the flux-functions A ± : 1R — > M d , with components (af)i=i t ... d, are assumed to be twice 
differentiable for definiteness. An initial data w(0,x) = wo(x) supplements this formulation, but 
obviously, some extra-condition, the coupling condition, must be prescribed at the interface H. 
For simplicity, we restrict ourselves in this introductory section to piecewise smooth solutions w 
with bounded left and right traces at the interface H: 

w(t,y±) = lim w{t,y±zv), y £ H. 

z— >Q+ 

Then, it sounds natural that the coupling condition we seek should relate these traces 

£(w(t,y-),w(t,y+))=0, t > 0, y £ H, (2.2) 
for some nonlinear mapping £ to be specified. The implicit function theorem is assumed to apply 



so as to recast (2.2) in the more tractable form 

w(t, y—) = c(w(t, y+)), t>0, y£H, 



(2.3) 



for some function c mapping 1R onto 1R. Assuming from now on c to be strictly monotone, we 
re-express the above coupling condition in terms of two nonlinear monotone functions 9- and 9+ 
with c = 9Z 1 ° 9 + : 



iw(t,y-))=9+{w(t,y+)), t>0,y£H. 



(2.4) 



Here and without loss of generality, and 9 + are assumed to be strictly increasing and to map 
1R onto -ZR, and their inverse functions are denoted by 7_ and 7+. On the basis of this pair of 
functions, we introduce the following useful change of unknown: 



u{t, x) 



9-(w(t,x)), t>0, x £ £>_, 
9+{w{t,x)), t>0, x £ V + , 



so that the coupling condition (2.4) resumes to 



u{t,y-) = u(t,y+), y£U. 



(2.5) 



(2.6) 



Observe that in the new unknown, (2.6) juste reads as a continuity condition for u. 



It is worth underlining that (2.6) defines the constant solutions of the coupling problem (2.1) 



(2.4), i.e. time independent functions w = w(x) which solve (2.1) and (2.6). Such functions clearly 



obey 

u(w(x)) = u*, x£M d \n, (2.7) 
for some real u* £ M. This observation actually just opens a path toward the mathematical study 



of perturbed solutions of the trivial solution (2.7). We refer the reader to the work [TT] devoted 
to the existence of self-similar coupled solutions for systems. 
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Observe that the coupling condition (2.4) plays the role of a pair of transient boundary con- 



ditions for the interface H. In other words, the coupling framework we address merely takes the 



form of two nonlinear hyperbolic IBVPs linked via the transient boundary condition (2.4). It 



becomes clear that the coupling condition (2.4) is actually expressed in a strong sense, since it is 



formulated without reference to the signature of the wave speeds at the interface H. It is never- 
theless well-known that the sign of the wave velocities at a boundary directly affects the boundary 



condition to be prescribed. Hence, the coupling condition (2.4) or its equivalent form (2.6) must 
be stated in a weak sense. 

We follow the approach for coupled problem in one space variables, originally developed in 
Godlevski, Raviart, and collaborators (cf. [371 [35] and 0H1IIS]). ^ n these papers, a weak form 
of the coupling condition (2.4) is formulated in terms of an admissible boundary set, proposed 



by Dubois and LeFloch [24] and based on the notion of Ricmann solutions. Such a notion here 
readily extends since the coupling condition expressed in (2.4) just links the traces of the coupled 



solution w in the normal direction v and thus essentially concerns the quasi-one dimensional form 



of (2.4) written for plane wave solution in the ^-direction. Thus it turns natural to consider the 



coupled problem in one space variable (up to some shift in the space variable z) 

d t w + d z A± (to) = 0, t > 0, ±z > 0, 



(2- 



where we have set A„ (w) — A (w) ■ v. In order to state the weak form of the boundary condition 
9-(w(t, y— )) = 6+(w(t,y+)), y G H, we first recall the Dubois-LeFloch framework for say the 
right IB VP: 



d t w + d z A+{w) = 0, t > 0, z > 0, 
w(t,0+) = b, t>0, 



(2.9) 
(2.10) 



for some prescribed real b. Following Dubois and LeFloch, a weak formulation of (2.10) is stated in 



terms of Riemann solutions associated with ( |2.9[ ), that is, W(-;wl,Wji) (for left- and right-hand 
states wl, wr): 

w(t, 0+) G 0+{b) = { W(0+; b,w),w G R}. (2.11) 



Observe that the analogous of (2.11) for the left IB VP built from A v would read 



;(t,0-) G 0~{b) = {W(0—,w,b),w G JR}. 



These considerations naturally yield us to the following coupled boundary conditions (2.4) at any 
point y G H and for t > 0: 



w(t,y+) e 0-(9^ o e^(w(t,y~))), 
w(t,y-) G Otiei 1 o 9 + (w(t,y+))), 



(2.12) 



This simple problem, based of two coupled equations at a given hyperplane, can be easily 
extended to more general interfaces resulting from a partition of lR d into two non-overlapping 
open sets T> + and X>_ such that 2?_ UV + = M d , separated by a smooth boundary 3D = 2?_ nl?+. 
Smoothness allows to define without ambiguity an unit normal vector v(y) for all y G 3D so that 
left and right traces at dD for piecewise smooth solutions of the coupled problem (2.4 1 may be 
defined as follows: 

w(t,y±) = lim w(t, y ± zv(y)), y € &D. 



The expected coupling condition just takes the weak form (2.12) 



Coupling technique based on an augmented PDE's system 

As already emphasized, an alternative coupling framework has been introduced by the authors in 
[IT] . Instead of dealing with two IBVPs coupled at a given interface via boundary conditions, our 
new approach treats the coupling problem as a single initial value problem, over the entire space 
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M d via an augmented PDE formulation. This strategy was introduced by the authors [ITJ [12l [13] 
for problems in one space variable. In order to encompass problems in several space variables, we 
perform hereafter a comprehensive derivation. 

The derivation starts from the characteristic functions of the two open sets 2?_ and 2?+, we 
denote by 

v- = Xx>- i v + = Xv+ ■ 



It heavily makes use of the change of unknown u introduced in (2.5), we rephrase as: 

\9 + (w(t,x)), ifv-(x) = 0, i.e. if u+(a;) = 1, 

Equipped with these notation, we recast the two distinct hyperbolic equations in T>± in terms of 
U: 

d 

-f' ± (u)d t u + ^i ± {u)at\^±{u))d Xi u = Q, t>0, x£V ± , 

8=1 

restricting ourselves to smooth solutions in a first stage. Recall that 7+ (respectively 7_) denotes 
the inverse function of 9 + (resp. #_). We further proceed by rewritting the above two equations 
in term of a single equation in x £ R d \ dD: 

d 

(v_~/'_(u) +v + "f' + (u))d t u + ^2 ( v -l'-( u ) a ih-(u)) + v + j' + (u)af('y + (u))Jd Xi u = 0, 

i=i 

At this stage, it must be noticed that the two characteristic functions u_ and v + in the above 
equation may be replaced by a single function say v, by setting for instance V-(x) = 1 — v(x) 
and v + (x) — v(x) for x £ M d \ dD with v — xv + - In the following, such a function v will be 
refered to as a color function. For the moment v is nothing but a step function taking values in 
{0, 1} but it is important to conceive v as a function taking values in the interval [0, 1] so that 
the value restores the equation set in T>_ while the value 1 restores the equation set in T> + . 
Intermediate values of v then may be thought as modeling a smooth shift from one problem to 
the other. Keeping this in mind we now recast the equations above in the form of an augmented 
PDE system with unknown u and v, for t > and x £ R d \ dD: 

((1 - v)i_{u) + vi + {u))d t u + ((1 - v)y_(u)WA-(^(u)) + t; 7 ;(7i)VA+(7 + (w))) • V x u = 0, 
d t v = 0. 

(2-13) 

We stress that the 1-dimensional form of these equations written for plane wave solutions in the 
direction v reads it > 0, x £ R d \ dD, or ±z > 0): 

((1 - «)tL(«) + wy' + (u))d t u + ((1 - v)j'_(u)VA-("f-(u)) ■ v + I)7 ;(«)V4 + ( 7+ («)) • v)d z u = 0, 
d t v = 0. 

(2.14) 

This system is easily seen to be hyperbolic if (and only if) the following quantity is not zero 

(1 - v)j'_(u)VA-{~/-(u)) ■ v + v-f' + (u)VA + {"f + {u)) -i/^O. (2.15) 

For such states, the standing wave associated with the additional unknown v can be seen to admit 
u as a Riemann invariant. In other words, as long as the non-degeneracy condition (2.15) is valid, 
u stays continuous at the jumps of the color function v, namely across the coupling boundary 
dD at which the value of v shifts from to 1. In other words and whenever (2.15) is valid, the 
coupling condition (2.6) is satisfied in the strong sense across the standing wave 

u(t,y-) = u(t,y+) ) y£dD. (2.16) 
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Violation of the condition (2.15) at a point of jump for v, namely at the interface dT>, expresses 



that waves from the left and right propagate with opposite sign at the interface; the first order 
system (2.14) is then only weakly hyperbolic. This is the resonance phenomena for which we refer 
the reader to, for instance, Goatin and LeFloch 126 and the references cited therein. As far as 



the coupling issue is concerned, the continuity condition (2.16) is no longer satisfied and the weak 



form (2.12) of the coupling condition must be addressed. Turning considering the augmented 



formulation (2.14), resonance phenomena has been studied in depth in [TT] in the scalar case 
thanks to a self-similar viscous perturbation. The Ricmann solutions for (2.14) defined in the 



limit of vanishing viscosity satisfy (2.12) when resonance takes place. To sum up, weak solutions 



of the augmented equations (2.14) and thus their multi-dimensional form (2.13) naturally encode 



the weak form of the coupling condition. 

We now generalize the rather special form of the augmented equation and adopt the general 
framework introduced by the authors in [TT] (which also applies to systems in one space variable). 
We thus introduce coupling functions Co : Mx [0, 1] — > JR and C, : Mx [0, 1] — > M with i G {1, . . . , d} 
satisfying the following consistency properties: 



\im v ^ Q C (u,v) = 7_(u), 
\im v ^ Ci(u,v) = a,r(j_(u)), 



lira,; 
lim,. 



nC (u,v) = 7+(u), 
n Ci(u,v) = a+(7 + (»), 



so as to consider in place of ( |2.13 1 the general augmented equations 

d 

d u Cp(u, v)d t u + d u Cj(u, v)d Xl u = 0, 

d t v = 0, 



i=i 



> 0, x e Sir, 



(2.17) 



(2.18) 



which equivalently recasts as: 

d d 
d t C (u,v) + ^2d Xz Ci(u,v) -y^d v Cj(u, v)d Xi v = 0, 



i=l 



t>0, x G ST 



d t v = o. 



In the following, the coupling functions Co and (Ci)i<j<d are smooth and 

C ,(Ci)i<i<d eC 2 {Mx [0,1]), 
and Cq, in addition, obeys 9„Co(u, v) > 0, u G M, v G [0, 1], which is a non-degeneracy condition 



for the time arrow in (2.18) 



The resonance phenomenon is the main difficulty in the coupling problematic and has made 
the matter of previous works especially in the one-dimensional case [TT1 IT2l ITS"] . In this one- 
dimensional setting, the analysis proves that if the resonance occurs for ( |2.18[ ) the self-similar weak 
solutions obtained via self-similar regularization satisfy the coupling relation (2.16). Nevertheless 
in the general case where resonance may appear, uniqueness then generaly fails for the initial value 
problem. 



The central interest of the augmented formulation (2.18) over more classical coupling ap 



proaches built from a collection of IBVPs stems from the fact it can be supplemented with a 
variety of regularizing mechanisms at the coupling interfaces. These regularization mechanisms 
are intended to handle the resonance phenomena which is likely to take place at the interfaces. A 
first regularization procedure relies on introduction of suitable viscous mechanisms. Such mech- 
anisms yield a non trivial internal structure to resonant interfaces which proves to be useful in 
the selection of discontinuous solutions. It turns that discontinuous solutions may not be unique 



for thin interfaces. The augmented formulation (2.18) actually allows for another regularization 



mechanism based on thick interfaces. The color function which is naturally discontinuous (for the 
description of thin interfaces) is regularized in the thick regime. Such a regularization technique 
has been analyzed in one space variable, and existence and uniqueness of a solution for the Cauchy 
problem has been established. In the next section, we show how to extend this regularization pro- 
cedure to several space variables. 



7 



Remark 2.1. An example of coupling functions satisfying the above conditions is 

C {u, v) = (1 — u)7_(ti) + U7 + (it), 

Cj(u, u) = (1 — u)a^(7_(u)) + ua+(7 + (w)), 1 < i < d. 

2.2 A framework for multi— component coupling problems 
Multi-component coupling of initial boundary value problems 

We are in a position to present the general coupling framework we intend to analyze in this paper. 
The proposed extension treats the coupling of (£ + 1), L > 1, distinct conservation laws in several 
space dimensions, with possible covering. The coupling modeling via augmented PDEs relies on 
a partition of the space lR d in a finite number of non-overlapping, non-empty and open sets 
(D%<i<l: 

L 

[JW = lR d . (2.19) 

1=0 

The set of boundaries B are given by 

B={JV k nV l . (2.20) 

An interface %m is by definition the part of the boundary of T>k which is only shared with T>i (see 
also Fig. [Tjfor an example with N = 2 and L = 3): 

Uu = (^H^)\ [J% (2-21) 

These interfaces Jiki are supposed to be smooth enough so that they admit an unit normal vector 
v ki{y), which is well-defined except at some "exceptional" points (like corners, etc.). We suppose 
the set of boundaries B to be of d-dimensional Lebesgue measure zero, and, more precisely, the 
remaining set B \ (Uk^iHki) nas only components of Hausdorff dimension less than or equal to 
(d — 2) (see for example the four points underlined in Figure [lj. 




Figure 1: Boundaries (in bold-face "Ho2, circle points being excluded) 



In each domain T>i , the unknown w is governed by a specific conservation law with flux- function 
A 1 = (<4)i<i< d :w6fi4 A l (w) G M d : 



dtw + ^2d Xi a[(w) =0, w(t,x)eJR, t > 0, leP, 



(2.22) 



Following the description introduced in the previous section, we start focusing the discussion 
on the definition of constant states ( 2.5 )-( 2.6 )-( 2.7 1 for the global problem set on the whole space 
M d . These solutions are recovered through a certain change of variable in each subdomain T>i, for 
Z = 0,...,L, 

u(t,x) =6t(w(t,x)), t>0,xeT>i, (2.23) 
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so that the stationnary solutions w(x) for the coupled problem (2.22) are the real constants u* in 
the u variable: 

u(w(x))=u*, xeM d \B. (2.24) 

The coupling functions 9i arc supposed to map increasingly 1R onto itself and we denote once 
again 7; the inverse functions: 

7l = flr X . 1 = 0,. ..,L. (2.25) 

Observe that a different outlook where the coupling functions would be associated to the interfaces 
T-Lki rather than to the domains themselves could only be local in space and therefore would not 
allow a matching of local constant solutions so as to define a global constant solution. Here we 
take advantage of the local formulation at each interface in terms of the traces of w, say w(t,y h ) 
and w(t,y) on the T>k~ and 2?/-side ofHki, respectively (relatively to its normal Vkl(y))'- 



°(w(t,y k ))=e l (w(t,y 1 )), t>0,yen 



kl- 



The following augmented PDE formulation is based on a vector-valued color function that 
merges the description of the coupling problem. In this multi-domain approach, this function is 
based on the set of characteristic functions of each domain: 



VL = XT> L 



so that the change of variable (2.23) may be also rewritten 

u(t, x) = 9t(w(t, x)), x e M d \B such that vi(x) = 1. 



(2.26) 



(2.27) 



Observe that since the (£ + 1) domains are a partition of the whole space M d , only L of the above 
characteristic functions are useful to complete the coupling description of the (L + 1) domains. 
Up to some relabeling we choose Vi,...,vl, so that vq is recovered thanks to 



vg(aO = l-X)ui(aO, xeM d \B. 



(2.28) 



1=1 



Multi-component coupling based on an augmented PDE's system 

In the following we make use of the vector-valued color function v = (i>i, . . . , Vl). At this stage, it 
takes values in the discrete set {0} U {ei} U . . . U {cl} where e/ stands for the l-th canonical vector 
of M L . This color function is intended to be regularized and to take values in the convex hull 
l£ = {v = (vi, . . . ,v L ) e R L /vi > 0, Ya=i ^ < !}• Tne problem \2.22) is then understood in 
the augmented form (with t > 0, x e M d ) 



d u C (u, v)d t u + ) j d u Cj(u, v)d Xi u = 0, 
»=i 

d t v = 0, 



(2.29) 



where the coupling functions Co and Ci are assumed to restore the formulation (2.22) in terms of 
u in each open set T>i, that is: 



lim„^ C (u,i;) =j (u), lim„_>. ei C (u, v) =7;(n), 
tirov-to Ci(u,v) = a°(7o(w)), hm„^ e; C t (u,v) = o|(7i(u)), 

The following smoothness and monotonicity assumptions are required 

d u C o {u,v)>0, ueM, «6lJ. 



1< i < d. 



(2.30) 

(2.31) 
(2.32) 



9 



This last property ensures the validity of the change of variable u i— > Co(u, v) for any fixed v, and 
the non-degenerate nature of the time- arrow in the augmented equations (2.29). 
In this context, the augmented system in the main unknown u reads 



d L 



d t C Q {u,v) + } j d x .C i {u, v) - } j d Vl Cj(u, v)d Xi vi = 0, 

d t v = 0. 



i=i 



i=i i=i 



(2.33) 



In the following, it will be useful to consider the same system written in the variable w — 
Cq(u,v) (denoted by w(u,v), and with inverse u(w, v) for each fixed v). Equipped with such a 



change of unknown, (2.33) becomes 

d 



d L 



d t w + ^2d x Ji(w,v) - l\{w, v)d x .vi = 0, 

d t v = 0, 



i=i i=i 



(2.34) 



where fi(w,v) — Ci(u(w , v) , v) and £\(w,v) — d Vl Ci\ u (u(w,v) : v) with i E {l,...,d} and I 6 
{!,...,£} (i.e. I — V V C). Hereafter and to shorten the notation, we write 



d t w + V • f(w, v) — £(w, v) : Vv = 0, 

d t v - 0, 



(2.35) 



with obvious notation. 



Entropy stability and well-posedness 

As already emphasized, in this work we propose a regularization mechanism based on thick in- 
terfaces that are modeled by any suitable regularized version of the discontinuous vector-valued 
color function v introduced in ( 2.26 )-( 2.28 ). For definiteness, we shall consider color functions v 
in W 2 >°°{1R + x R d ,M\). Obviously, it suffices to choose the initial data v in W 2 '°° (M d , Bjj s o 
as to inherit from the required smoothness in the v solution of the augmented equations (2.35). 
In turn and arguing about this smoothness property, the equations under consideration reduce to 
an inhomogeneous scalar equation in w: 

d t w + V ■ f(w, v(x)) = £(w, v(x)) : Vd(i), (2.36) 

where the right-hand side just plays the role of a classical source term; namely this term does not 
contribute to the definition of the possible discontinuities of w. At a point of jump, (2.36) just 
resumes to the classical Rankine-Hugoniot condition 



a(w+ -w~)+J2 (fi(w+,v) - fi(w-,v)) = 0. 



(2.37) 



A selection criterion of the admissible weak solutions w is of course needed, and we recast the 
balance law (2.36) in the main variable u: 

d t C (u,v) 



d 

£ 

i=l 



d u Ci(u,v)d Xz u = 



(2.38) 



for all smooth solutions. For such solutions, additional equations are deduced and based on any 
(strictly) convex function w n- U(w), by multiplying (2.38) by U'(Cq(u, v)), 

d 



d t U(C (u,v)) +^2d u Q l (u,v)d Xi u = 0, 



(2.39) 



i=l 
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where 



Qi(u,v) 



U'(Co{6,v))deCi(9,v)d6, l<i<d. 



We thus get from (2.39) the equivalent form for smooth solutions u: 

d d L 



d t U(C {u,v)) + ^2d Xi Qi(u,v) = ^2^2d Vl Qi(u,v)d x 



vi- 



(2.40) 



(2.41) 



Observe that the above right-hand side is nothing but a classical source term since we again 
emphasize that the color function v is smooth. As a consequence, the weak form of (2.411 for 
discontinuous solutions u reads: 



(u,v) < y^y^<%, Qi{u,v)d x 



vi, 



(2.42) 



which naturally plays the role of an (inhomogenous) entropy inequality for selecting the relevant 
weak solutions. Hereafter, we shall make use of the inequalities (2.42) for all convex entropy U. 
These will be alternatively invoked (essentially when the color function is locally constant) in the 
w variable: 



d t U{w) + dx^iiw, «) - S C '^ w ' v ^ Vl - °- 



i=i i=i 



with 



(2.43) 



Fi(w,v) = Qi(u(w,v),v), Ci(w,v) = d v Qi\ u (u(w,v),v), 



Ki< d. 



To shorten the notation, equation (2.43) are written as 

d t U{w) + V ■ T(w, v) - C{w, v) : Vt> < 0. 



(2.44) 



(2.45) 



The inhomogeneous scalar conservation law (2.36) supplemented with all the entropy inequalities 
(2.43) naturally falls within Kruzkov 's theory of entropy solutions, since the color function v 
belongs to W 2 - °(]R d ,M^). Therefore, Kruzkov's uniqueness theorem for scalar conservation law 
with smooth inhomogeneities applies and asserts the uniqueness of the entropy weak solution of 

the Cauchy problem ( |2.36[ )-( |2.43[ ) with initial data w € L 1 (5? d ) n L°°(lR d ). 

Hereafter, we shall prove existence and uniqueness of a solution to the coupled problem ( 2.36 )- 



(2.43) thanks to a multidimensional well-balanced finite volume method formulated on general 
triangulations. Here, the well-balanced property means that the solutions in the u variable is kept 
constant in time and space as soon as the initial data uq is chosen constant whatever the definition 
of the (smoothly varying in space) color function v. This well-balanced property is obviously a 
constancy property of primary importance. 



3 A well-balanced finite volume scheme for coupling prob- 
lems 

3.1 Terminology and assumptions 

Before stating our main result, we introduce some notation and motivate our formulation of the 
finite volume method under consideration. To meet the well-balancing property, the finite volume 
framework we develop uses two families of triangulations. The first triangulation, denoted by 7h, 
is made of general polyhedra and will be referred to as the primal mesh. Then a closely related 
triangulation is of concern, the dual mesh 7ft*, whose polyhedra are derived from the edges of the 
primal one. As we shall see, dual meshes may not uniquely defined from Th and it will turn that 
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a given choice essentially affects the closed-form of expression of the CFL restriction in the (time 
explicit) finite volume method. 

Equipped with these primal and dual meshes, approximate solutions Uh and Vh of the Cauchy 



problem (2.33) with initial data (ho,i>o), are sought as piecewise constant functions. In constrast 
with the usual approach, constant values for Uh and Vh will not be co-localized: Uh (and Vh, 
respectively) will assume constant values in each polyhedron of the primal mesh (and the dual 
mesh, resp.). 

To facilitate the derivation of the proposed well-balanced scheme, we shall take advantage of 
the regularity of the color function v, which provides some room for the specific definition of the 
discrete approximation Vh- it may range from a local averaged form to a point- wise evaluation. 
Without real loss of generality, we use an average value of v along each edge of the primal mesh. 
This choice allows to bypass the definition of the dual mesh from the edges of the primal one: a 
convex sequence of reals, in turn, provide sufficient information on the dual mesh. On the ground 
of this observation, we shall give a first brief but sustained mathematical presentation of the finite 
volume method under consideration. We shall then be in a position to state the main result of 
this paper. At last, we shall close this section with a comprehensive construction of the proposed 
finite volume approximation when deriving dual meshes from the primal one. 

The primal mesh, Th, is a general (locally finite) triangulation of M d made of non-overlapping, 
non-empty, and open polyhedra : \JKeT h K — JR d - We assume that for every pair of distinct 
polyhedra K, K' E Th the set K n K ' is either an edge e of both K and K' or a set with Haussdorf 
dimension less than or equal to d — 2. The set of edges of a polyhedron K is denoted by dK; 
and for each e € dK, vk.c € JR d represents the outward unit normal vector to the edge e (see 
Figure [2]). The volume of K and the (d — l)-measure of e are denoted \K\ and |e|, respectively. 
Given an edge e in K, K e denotes the unique polyhedron in Th that shares the same edge e with 
K. We set h = sup Ke7 - h Kk, where fix is the exterior perimeter of the polyhedron K, and assume 
that the triangulation Th satisfies the following non degeneracy condition 

h-KPK ^ r , fn -J \ 

sup ~T^T - C ( 3>1 ) 

for some constant C > 0. Here, px denotes the perimeter of K defined by px = Seeax M- 

It is unnecessary, at this stage, to provide a comprehensive derivation of the dual mesh Th* that 
one could define from the edges e in the primal mesh Th- Recall that, by design, a dual mesh is made 
of non-overlapping, non-empty, and open polyhedra denoted by K*(e) with U e( =j- h K*(e) = IR d . 
By construction, both sets K*(e) n K and K*(e) n K e are non-empty for all pair (K,K e ) of 
adjacent polyhedra parametrized by the edges e in Th- Note that the set K*(e) n K is a subcell of 
K. Then, the only information about Th* that is required in this section is a given convex sequence 
of reals prescribed in each polyhedron K in Th', we denote by {ctK,e}{ e .eedK}, that satisfies (for 
any K in Th) 

< a K ,e < 1 (e € dK), a K,» = L ( 3 - 2 ) 

eedK 

We will see later that the coefficient Q.k,b is nothing but the ratio of the volume of K*(e) D K to 
the volume of K, where K*(e) stands for the dual polyhedron of K attached to any edge e in dK: 

\K*(e)nK\ 

a K ,e = r^n , e G dK. (3.3) 

At last, the time increment, denoted by r, is assumed to satisfy f < C and the primal mesh to 
be constrained by 

< M < C2 (3.4) 
n 

for some constants C,C\,C2 > 0. Whereas the latter is probably not an optimal condition, it 
sufficies to ensure the non degeneracy of the mesh: all one-dimensional characteristic lengths are 
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of order h. A key property for the forthcoming CFL condition, is that under these assumptions 
the area \K*(e) n K\ is not smaller than 0(h 2 ): there exists a positive constant c such that 

ch 2 < \K*(e)n K\. (3.5) 

We use the notation t n = tit. As already underlined, we will seek at each time level t n 
approximate solutions Uh and Uft of the Cauchy problem (2.33) with initial data (uq,vq), under 
the form of piecewise constant functions with: 

u h (x,t n ) = u n K , xeK, KeT h , 

(3.6) 

Vh(x,t n ) = Vh(x) = v e , x€K*(e), e e 7ft. 



Here and since the solution v in the Cauchy problem (2.33) does not depend on time, it seems 
natural to set Vh(x, t n ) = v(x) — Vf l °(x) € 1R L for all time level t n , for some discrete approximation 
Uft of the smooth function v . We introduce 

v h (x)=v e = ~7~~7 v a (y)dy, x€K*(e), e € 7ft, (3.7) 

M Je 

while the discrete version of the possibly discontinuous initial data uq is chosen according to the 
usual full averaging procedure over each polyhedron K: 

u h {x)=u K = —J^u {y)dy, x e K, Ke%. (3.8) 

Remark 3.1. Since vo is regular, any other consistent definition for the constant value v e in K*(e) 
would have been relevant. The interest in the particular choice {3.7) stems from the following Green 
formula, valid for each polygonal domain K : 

X v el v K>e \e\ = / V ■ (vi(x) X)dx = X / Vvi(x)dx, 



K 



where X denotes any fixed vector in lR d and v ei (and Vi, respective ly) t he l-th component of the 
vector v e € 1R L (and v, resp.). Hence the proposed average value in (3.1) comes with the identity: 
J K Vvi(x)dx — v el VK,e\&\- In a tens orial notation, we thus get J K Vv(x)dx = v e ®VK,e\&\- 

eedK eedK 

The evolution in time of the discrete solution uh will rely on a family of numerical Qux- functions, 
associated with each edge e of any polyhedron K in 7ft. Besides other properties, these numerical 
flux functions must meet some consistency property with the exact equation for governing u in 
(2.35), namely: 

d t w(u,v)+V ■ f(w(u,v),v)-£(w(u,v),v) : Vu = 0, x e K, te(t n ,t n+1 ). (3.9) 

Observe that in the neighborhood K*(e) of each edge e, where reduces to a constant value v e , 
the above equation boils down to the scalar equation in the unknown w = w(u, v e ): 

d t w + V • f(w,v e ) = 0, xeK*(e)f)K, te(t n ,t n+1 ). (3.10) 

This in turn leads us to define the required numerical flux function at each edge e in 7ft as a 
locally Lipschitz continuous two-point flux-function g e ,x{-i A Ve) ■ Si X Ml ]R that satisfies the 
consistency property: 

g e ,K(w, w; v e ) = f(w, v e ) ■ v K<e , (3.11) 

the conservation property: 

g e ,K{w, w e ;v e ) = -g etKa (w e ,w;v e ), (3-12) 
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for all reals w and w e , and the monotonicity property 



dg(w,w e ;v e ) 
dw 



>0, 



dg(w,w e ;v e 
<9w. 



< 0. 



(3.13) 



In addition, we assume that the numerical flux depend (locally) Lipschitz continuously in the 
variable v e . 

Standard 3-point monotone schemes in the scalar framework obey ( 3.11| )H3.13) and that the 
main results in this paper are easily extended to all E-schemes (Oshcr [40 ). For clarity, the 
dependence in the parameter v e appears explicitly in the numerical flux-function g e ,K{-, ■', v e ). 

Remark 3.2. Since the function $(.,.;.) : M 3 — > 1R is locally Lipschitz continuous in its three 
arguments, for all compact K. C M 3 , there exists some positive constant Ck: such that for all triple 
(w^,We ,vi ) and [w^\w^\ve) infC, the following estimate holds true: 



\g e , K {w^\w^;v^) 



9e,K 



{wP>,wP;vP)\ < C K (\w^ - W W\ + \ w f) - W V)\ + \ v f) - U (D|). 



3.2 Definition of the well-balanced scheme 



We are now in a position to define the finite volume approximation of (3.9). Assuming that the 
approximate solution Uh{-,t n ) is known at time t n , we determine the evolution up to the next time 
level t n+1 as follows: 

Subcell reconstruction. At each time t n in each polyhedron K of Th, we consider any edge e G OK 
and introduce the subcell state 



WK,e = Co(uK,Ve), e € OK, 

as well the following average over all edges of K 



eedK 



OiK,eW K , 



(3.14) 



(3.15) 



Evolution in time. In order to the discrete solution Uh at time i n + , we define (in each polyhedron 
K) to be the unique solution of 



^2 UK.e C (u^ +1 ,w e 

eedK 



- w n+1 



(3.16) 



where the state w^ 1 " 1 is given by the finite volume scheme 



m n+1 — in 71 — 
W K — W K 



•K 



\K\ 



£ 9e,K{ W K,e> W K e , e ; V e)\ e \ + |4 E f( W K,e^e) ' V K 



1 eedK 



(3.17) 



1 eedK 



This completes the description of our numerical method. The proposed finite volume method 
is explicit in time and, for the sake of stability, we need to impose a CFL (Courant, Friedrichs, 
Lewy) condition which reads, for all polyhedra K in Th and edges e € dK, 



— sup 

\Jy I CtK,e ue[m,M] 



df(w(u,v e ),v e ) 
dw 



< 1, 



(3.18) 



where m — inf uo(x) and M — sup uq(x). 
xeR d x£lR d 



Due to the dimensional hypothesis (3.2 )-(3.4|-(3.5 1 the ratio \K\ax te /\e\ satisfies 



a^ e \K*(e)nK\ c 
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so that the CFL condition can not imply the degeneracy of the time step r, that decreases at most 
as O(h). We will see in Section [I] how to build suitable primal and dual meshes 



Several comments are in order. First observe that the constitutive assumptions ( 2.31 |— ( 2.32 1 



on the coupling function Co(., .) immediately yields existence and uniqueness of a solution to the 
nonlinear equation (3.141 so that the finite volume method ( 3.14 )-( 3.17) is well defined. The 



formulas (3.14) and (3.16) obviously express the same identity at the times t n and i n+1 , and are 



redundant: the finite volume method essentially reduces to (3.16H3.17). As they stand, they 
nevertheless ease the description of the method. 



Next, it is worth observing that the consistency condition ( |3.11 1 allows in (3.17) to recast 



the flux balance > f{w 1 j { el v e )vfc,e\e\ as ^ ^ 9e,K( w K ,ei w K,e> v e)\e\. Here we stress that at 



v e ) and its counterpart 



eedK eedK 

each edge e in dK, both the numerical flux- function ge.ir^Ki 
f( w K.e> v e) ' v K,e are evaluated thanks to the subcell values w K e (3.14) and not to their averaged 
for m w 1 ^ in ( 3.15[ ). The motivation is twofold. In a first hand, the two flux balances involved 
in (3.17), namely 5e,if| e l and f{w r j <e ,v e ) ■ VK.e\e\, make the proposed formula to be a 



eedK 



eedK 



consistent finite volume approximation of the exact equation (3.9 1 for governing u: namely, the first 



one will be seen hereafter to be consistent with V • f(w, v) while the second one actually provides 
a consistent approximation of the source term £(w, v) : Vw. In a second hand, the discretization 
of the source term is seen to be well-balanced. 



Proposition 3.3 (Well-balanced property). When the initial data uo for \3. 9\l is a constant 
function uq(x) = u*{x € lR d ), then, for any choice of the color function v in (3.9), the discrete 
solution Uh of ( 3.1J^ -(3.11) is also constant, with 



u h (x,t n ) = u (x) = 



(3.19) 



for all time level t n . 



natural equilibria of (3.9) 



In other words, the finite volume method (3.14)-(3.17l is well-balanced with respect to all the 



u* for all x in 1R so that at the 
for any edge e of an arbitrary 



J K c ,e 



Proof. The discrete initial data (3.8) clearly reads uf^x) 
first subcell reconstruction step, we get uP K e — Cq(u*,v c ) = 
polyhedron K in 7ft. Consequently, the numerical flux g e ,K{w^ e ,ffl^ e ;w e ) at any edge e boils 
down to f{w° K e , v e ) -VK.e in view of the consistency condition (3.11 ). Namely the two flux balances 
in the updating formula (|3.17|) cancel out and we end up with w\r = w 



K 



^ aK, e C (u*,v e 



eedK 



thanks to the definition (3.15). Arguing about uniqueness, we thus get when solving ( |3. 16 ) v} K = u* 
for any polyhedron K of 7^: namely Uf l (x, t 1 ) = u* for all x in M d . An immediate recursion extends 
the result to the subsequent time levels. □ 

To conclude this paragraph, it is worth illustrating that the last flux-balance entering the 



finite volume approximation (3.17) actually provides a consistent approximation of the source 



term £(w, v) : Vw. For the sake of simplicity, we temporarily adopt (cf. Remark 2.1 ) 



C (u,v) = (1 - u)7-(«) + vj+(u), 

C t {u,v) = (1 - v)a~(7_(M)) + ua + (7+(u)), 1 < i < d, 
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so that f(w,v) and £(w,v) in ( 3.9 ) read f(w(u,v),v) = (1 — v)A (j_(u)) + vA + (j + (u)), and 
£(w(u,v),v) = (^A + (7_|_(u)) — A~(7_(it))J. It can be then readily computed: 

E f( W K,e, V e) ■ VK,e\ e \ 
eedK 

= E ((l-Ve)A-( 7 .(u n K ))+v e A+( 1+ (ul))) -v K M 

eedK 

= (V( 7+ «)) - A-( 7 _(«£))) • J] « e |e|i/x, e + A-(7_(«^)) ■ ( £ |e|^„ 

which is nothing el se a consistent discretization of £(w(u,v)) : V«, in view of the representation 
formula in Remark 



3.1 



(for Vu) and the identity | e| Kpc,e = 0. 



e<£dK 



These straightforward calculations allows to bridge the finite volume formula ( 3.17[ ) to the 



governing equation (3.9 1 for w(u,v), expressed over K, namely where Vh does achieve distinct 



values. The gap in between (|3.9[) and its reduced version (3.10) (i.e. with x € K*(e) D K) will be 



definitely closed when revisiting the finite volume approximation ( 3.14 1 — < 3.17 ) with primal-dual 
meshes (in Section [4]). 

3.3 Main convergence result 

We are now in a position to state the main result of this paper. 

Theorem 3.4 (Well-balanced finit e volum e method for multi-dimensional coupling problems). 
Consider the Cauchy problem fjj^-fjjjl) with initial data u a e L°°(R d ) and v e W 2 ' ao (lR d ) 



under the constitutive assumptions {2.31) 



2.32). Let Uh be the sequence of approximate solutions 



defined by the finite volume method \3.1)-(3.8) and (3.14-)~(3.lD with numerical flux-functions 
satisfying the conditions \3.11 )-(3.13). Then under the CFL restriction (3.18), the sequence 
is uniformly bounded in L°°(]R+ x ]R d ) and converges (when h 0) in the L\ oc norm strongly 
(1 < p < oo) to the unique entropy solution u to the problem (2.33)~ \2.J^2 ): namely for all time 
T > and for all compact K. in lR d 



lim \\u - Uh\\LP((o,T)xic) 



The rest of this paper is devoted to a proof of this theorem. 



4 Finite volume approximations on primal-dual meshes 
4.1 A convex combination 

One of our objectives in this section is explaining how the coefficients ctK,e should be determined. 
Arguing about the formula-definitions (3.14 H 3.15 ) at time t n and the consistency condition ( 3.11 ), 
we obtain the following statement. 

Lemma 4.1 (Edge values and convex combination). For any polyhedron K of Th and edge e in 
dK , let us define the following subcell states: 



ra+1,- 
' J K,e 



a K ,e \K\ 



.K 



/ 77 77 \ /77 77 \ 

( W K,ei W K^ V e) ~ 9e,K(W K W R ] V e ) 



(4.1) 



Then in (3.17) are recovered by the following averaging procedure: 

n+l 71+1.- 
W K = 2^ a K^ W K,e ■ 
eedK 



(4.2) 
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Observe that the finite volume formula (4.1 ) for w 



n+l,- 
,J K,e 



is nothing but a consistent approxima- 



, n+l, 



tion of the one dimensional conservation law: d t w + V • f{w, v e ) — 0. The reason for calling u>^- + e 
a subcell state will be explained in this paragraph and is at the core of the re-interpretation of the 



finite volume formula (3.171 with primal-dual meshes. 



To further proceed, let us underline that the identity (4.2) just expresses that w 1 ^ 1 actually 



is a convex decomposition of the subcell states w 



When understood in their quasi-one 
dimensional form (4.1 1, the latter can be recognized as extensions to the present inhomogenous 



n+l, 



setting of partial states entering similar convex decompositions that have proved well suited in 
the analysis of homogeneous multidimensional finite volume methods [T51 [T5] . Indeed, the interest 
in such a convex decomposition primary stems from the fact that many of the basic stability 



properties satisfied by the scheme (4.1) in one space variable are right away inherited in several 
space variables thanks to convexity under some CFL restriction. Observe that the relevant CFL 



condition for (4.1) reads 



\K\ a K . 



J K e ,e' 



Ve) - g e ,K(w 



K.e' w K,e 



>;0 



J K,, 



< 1, 



(4.3) 



and hence the CFL restriction (3.18) 



At last and arguin g abo ut the definition ( |4.1| , the subcell reconstruction step (3.15) at time 
i n+1 and the formula ( 3.16[ ), we deduce the (seemingly trivial) identities 



eedK 



En+l. 
a K,eW Ke 

eedK 



(4.4) 



In other words, all the steps involved in the method are locally conservative: this natural property 
will play a central role in the forthcoming analysis. 



4.2 A reformulation of the scheme 

The derivation of a dual mesh Th* from the edge of the primal one Th may be performed as follow. 
For any (open) polyhedron K, the idea is to pick an internal node xk in K which choice is left 
arbitrary at this stage. Such a procedure is given below a systematic definition independent of the 
mesh refinement h. Equipped with the node xk, we define for any edge e in K the convex hull 
of e and xk- The interior of this convex hull, we denote by £(xK,e), yields a non-empty open 
polyhedron made of (d+ 1) edges. Observe that the following properties are met by construction: 
for any pair of distinct edges e, e' in dK with K an arbitrary polyhedron in 7/j 

S(x K ,e)DK = £(x K ,e), £(x K ,e) r\S(x K ,e') = 0, (4.5) 
while £{xk, e) = K, Then, the required definition of the polyhedron K*(e) of the dual mesh 

eedK 

Th* , attached to any edge e in 7h with adjacent polyhedron K and K e , follows from 

K*(e)=£(x K ,e)U£(x Ke ,e). (4.6) 

We refer the reader to Figure [2] for an illustration. 

The constructive procedure for defining the internal node xk independently of h relies on the 
set of vertices $ of the polyhedron K, together with a convex sequence of reals {Pk^w^^k} 
satisfying: 

< K j < 1, ^J(; Yl Pk,* = L 

The required internal node xk in K is then defined by its coordinates in M d : xk = X^git Pk,S %■&■> 
where x® stands for the coordinates of the vertex <&. This construction ensures the correct behavior 
of the primal and dual meshes with the definition of the UK.e and with the previous non-degeneracy 
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Figure 2: Primal and dual meshes, edges and vertices. 



assumptions (3.3)-(3.5), the CFL condition (3.181 is then only modified according to the choice of 
the function v and its discrete representation. 

To further proceed in the comprehensive derivation of the finite volume framework, some 
additional notation is in order. For any K in Th and e in dK, an edge of a dual polyhedron 
K*(e) € Th* or of the subcell K*(e) fl K of K will be indifferently denoted by e*. Observe that 
with little abuse in the notation, an edge e of some cell K of the primal mesh Th is also a dual 
edge of the subcell K*(e) n K: see indeed Figure [2J At last VK*(e),e* € ^ d stands for the outward 
unit vector normal to the edge e*. 

Equipped with these not atio n, we are in a position to re-interpret the quasi-one dimensional 
state w^- + e 1,_ introduced in (4.1) in term of a state in the subcell K*(e) D K of K, thanks to the 



following simple but key identity: 

\e*\"K*(e),e* = °. 



E 



i.e. \e\v K , e 



2^ l e \ V K*(e), 

e*<£K*(e)nK, e'/e 



It is then straightforward to recast w^ e ' according to 

n+l,- ». T 



T 

E 

e*GA'*(e)ri-R", e*^e 



9e,K{wK^w^ ce ]v e )\e\ + 
9e,K[w n K ^w n K 



OL K ,e\K 

ON 



f(w K J ■ vk,> 



(4.7) 



f(w n Ke ,v e ) ■ u K *( B)te *\e*\), 



where we have used the interpretation (3.3| of ax, e - Introducing the numerical flux formula: 



9e*,K*(e) 



9e,K(WK,ei W K c ,e: v e), ^ e* = e; 

J( w K,ei v e) ■ ^K*( e ),e* , otherwise, 



(4.8) 



lo^g ' thus reads 



w 



K.i 



W 



K.i 



\K*(e)r\K\ 



E 9e*,K*(e)\e*\. 



(4.9) 



e* £K* (e)llK 



We can clarify the origin of the definition g e *,K*(e) = f( w Ke' v e) ' v K*(e),e* f° r edges e* distinct 
from e. For such an edge e*, it is worth introducing the adjacent subcell K*(e') f)K to if*(e) n K 
in X: i.e. with e' in Si-T such that K*(e') fl if*(e) = e*. Note that e* is of course distinct 
from e' . We then successively rewrite the left- and right-hand numerical flux at e*, say g e * t K*(e) 
(respectively g e * i K*{e'))i as follows: 



f{w{u r i i ,v e ),v e ) ■ VK*{e),e*, respectively : - f{w{u\,v e '),v e ') ■ v K *(e), 
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since by definition (3.14) uf^ e = w(u'] (: , v e ) and wk,c> — w(u^-,w e ') and, equivalently, 

(f(w(u,v),v) ■ V K *(e),e*)(u(0-)), 

respectively : - (f(w(u,v),v) ■ v K *{e),e* ) (^(0+)), 



(4.10) 



where u>(0 T ) stands for the left and right traces at £ = of the self-similar function u> : £ € 5?^ 
0(0X0) eMx M L given by 



w(0 



(t*5-,W e ), £<0, 

(«5-,« e '), C>o. 



(4.11) 



From Section [2j recall that the Riemann solution of 

<9 f W + d x (f(w(u,v),v) ■ V K *{e),e*) («, «) ~ d v (f(w(u, v) , v) ■ V K *(e),e 



Vv = 0, 



(4.12) 



(with initial data ((tt^-,u e ), a; < 0, {u'j i ,v e l ), x > 0)) consists in a standing wave separating 
(UfciVe) from (u^,v e '), and thus coincides with w(£) in (4.11). It is therefore clear that the 
flux-functions in (4.10) actually results from the Godunov method applied to the augmented 
system (4.12) at the edge e*. In other terms, the finite volume formula (4.8 )— (4.9 1 in each subcell 
K*(e) D K may be understood as an approximation of the balance law for governing u in (3.9): 

d t w{u,v) + V - f{w{u,v),v)-t{w{u,v),v) : Vu = 0, x e K, te(f\t n+1 ). 



This interpretation closes the gap in between the governing equation (3.9) for u and its reduced 
form (3.10) expressed in w: 

d t w + V ■ f{w,v e ) = 0, xeK*(e)DK, te(t n ,t n+1 ). 



4.3 Sup-norm estimates 

Throughout the upcoming sections, the assumptions of Theorem |3.4| are tacitly assumed to be 
valid. Their formulations are thus skipped over in any forthcoming statements. The main result 
of this section ensures that the sequence of approximate solutions Uh stays uniformly bounded in 
L°°{M + x ]R d ) as a consequence of the following result. 

Proposition 4.2 (Maximum principle). The finite volume method satisfies the following inequal- 
ities (in the variable u): 



mm | u K , mm u 



eedK 



K, 



< ut +1 < max 



e£dK 



'A', 



(4.13) 



in each polyhedron K in Th and at all time level t n 



Since vq € W 2,oc immediately implies a sup-norm estimate for Vh given by (3.7), we easily 
deduce, from the maximum principle (4.13), an additional uniform sup-norm estimate but for 
Wft = Co(uh,Vh) arguing about the regularity properties (2.31) of Co: 

\\w h \\ L ^ (M+xMd) <0(l). (4.14) 



Besides the monotonicity assumption (3.13) met by the numerical flux functions, we stress that 
the preservation of conservativity (4.4) in the subcell reconstruction procedure plays a central role 
in the validity of the reported maximum principle, as highlighted in the proof. The latter will be 
carried out using a recursion procedure based on subsequent partitions of the set of edges e in 
K. To fix the notation and up to some relabeling, {ei, . . . , ej K } represents the full set of edges 
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e G <9/T so that here the index Jk is given by #{e, e € dK}. Subsets of the form {ei, ej}, with 
increasing index J G {1, Jk}, will be of concern as follows. Being given J with 1 < J < K, let 
us attach to the subset {ei, ej} the solution u^rj^ ej } of the following nonlinear equation: 



H a K,e j Co( U K^{e x ,..., e jr V ej) 
l<j<J 



1<j<J 



4.1 



(4.15) 

Again, the constitutive assump- 
tions ( |2.31[ )-( |2.32| ensure existence and uniqueness of a solution to (4.15). 

Arguing about the conservation property (4.4) satisfied at the subcell reconstruction step, it 
is worth observing that uV&}~ , can be identified with the final state at time t n+1 



where the subcell states w K e . are defined in (|4.1|, Lemma 



in the finite volume approximation (3.14)-(3.17). Therefore, the recursion under consideration 
naturally ends up as soon as the index J reaches the value Jk- In order to initiate the recursion 
and propagate it, we need the following statement concerned with the values w^r e ~i j 1 < J < Jk, 

solutions of Co(«xT{ej}> We j) = W K^ej '• 

Lemma 4.3 (Local maximum principle). The maximum principle holds true at any edge ej in 
dK: 

mm(u n K , u n Kc j ) < } < max« , u n Kej ), 1<J<J K . 

Then the maximum principle "propagates" to sets {ei, ...,ej}, as follows. 



Lemma 4.4. The solution e ^ to (4-15) with J G {1, Jk}> obeys the following maxi- 

mum principle: 

min (u K , mmjtfk^ )) < u 7 ^'^ ej} < max (u K , mffl(u^ )J . 



The proposed lower and upper bounds for u 



n+l- 



K,{ ei ,...,ej K }> 



the estimate in the lemma with 



J = Jk, just reads the expected local maximum principle (4.13) for since again 

coincides with u n ^}~ , by construction. 

Proof of Lemma \4-3\ To alleviate the notation we skip the index J and first point out an estimate 
valid under the CFL restriction (3.18) for any edge e in dK: 



(4.16) 



as a well-known consequence of the monotonicity assumptions (3.13) satisfied by the numerical 
flux function g e ,K(-, ■', v e )- We then recall that the subcell reconstruction step (3.14) builds w K = 
Co(u K ,v e ) while the identity w^- + e 1_ = Co(u"^r~l,v e ) holds from our definition. We can thus 
recast (4.16) as: mm(C (u K ,v e ),Co{u Ke ,v e )) < C (u^~y,v e ) < max(Co(u K ,v e ),Co(u Ke ,v e )), 
from which we immediately deduce the desired estimate, namely 

min(M^-, u K J < u K ^~y < max(u Kl u K J, e G dK 



since the function Cq is by assumption (2.32) strictly increasing in its first argument. □ 

Proof of Lemma \4--4\ The desired lower-upper bounds with J — I are stated in Lemma |4.3| Then, 
assuming the validity of the maximum principle at rank </, 1 < J < Jk , this one is proved to hold 
at the rank (J + 1) starting from (4.15): 

1<J<(J+1) 

„n+l— 



a K,e^K + e)~ 
l<j<J 



OtK, 



e(j+l) W K,e (J+1) ' 



a K,e 3 C Q {u n +l ei ejV V e] ) + aK.e (J+1) C Q (u n K +^ j+i} ,V e{J+1) ] 

l<j<J 
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We recast the above identity as follows: 

X! a ^ C 0«t{eT,..,e (J+1) }'^ 3 )- Yl a *.^( tt £{e7,...,e J }> t, «i) 

l<j<J 



l<j<J 

= - at K . 



To condense the notation, we introduce the two functions u n- 9j(u) = aK, ej Co(u, v ej ) and 

i<j<j 

u !-> ^(J+i) ( U ) = a -ftT,e (J+1) Co(u,?; e(7+1) ) so as to deduce: 

...,e ( , +1)} ) ~ ^« { e" (^ + D(«5rfc,....^ +1)} ) - ^ + l)(«Stg7 +1) )) < o, 

since by assumption (3.2) % ie(; . > 0. But the monotonicity hypothesis (2.32) on Cq together 
with again assumption"! 3.2) imply that both functions u <— > ^>j(u) and u t-> strictly 
increase with u so that the above inequality yields 



Lemma 
pleted. 



4.3 



\ if,{ei,...,ej}' ft,ep +I) / - if ,{ei ,. . .,e ( J+1) } — V. if ,{ei ,e j} ' if ,e(j+i) > 

implies min(u^ r ,u^- ) < u^y"^ < max(u^-,u^- ), and the proof is com- 



□ 



5 Entropy inequalities 

5.1 Preliminaries 

Proposition |4.2| asserts sup-norm boundedness for the sequence which in the absence of an 
a priori strong compactness argument, leads us to study the structure of the Young measure /x 
associated with {uh}h>o- Recall that such a Young measure represents all the composite weak-star 
limits a(uh) of Uh with continuous functions a G C°(1R), namely for all continuous functions in a 
single variable 

a(u h ) — ^ <fj,,a>= / a(X)dfj,(X), 

JlR. 

weakly-star in L°° . We propose to establish that the measure /i under consideration reduces to a 
Dirac measure, and hence to prove the strong convergence of u^, invoking DiPerna's uniqueness 
theorem [22] for entropy measure-valued solutions. 

In this section we derive the required discrete entropy inequalities together with the a priori 
estimates that are needed to handle the passage to the limit in the sense of measure valued 
solutions. In this respect, the main issue is to assess the relevance of the Young measure /x in 
such a limit. Indeed, discrete entropy inequalities generically involve numerical flux functions, 
that are continuous functions but of (at least) two arguments: the sequence Uh{.) itself and its 
shift AhUh = Uh{- + h). Nonlinear superposition of possible discrete oscillations in ut and its 
shift AhUh may prevent the usual Young measure \i to represent the composite weak-star limit of 
G(uh, A/jU/j). Counterexamples have been constructed in Coquel and LeFloch [19] . Some weak 
control over possible discrete oscillations is therefore mandatory in order to justify the applicability 
of /i in the limiting form of discrete entropy inequalities. 

The requisite weak estimate corresponds to some estimate of the discrete entropy dissipation 
rate in the finite volume approximation. The derivation of several specific estimates with dis- 
tinctive features have been the matter of a large literature following Coquel and LeFloch [18.. 
(The reader is referred to the introduction where several subsequent contributions were quoted.) 
The estimates we derive now generalize the ones in Cockburn, Coquel, and LeFloch [TB]. The 
entropy dissipation estimate does not allow actually to pass to the weak limit in arbitrary numer- 
ical entropy-flux functions, but nevertheless turns out to be sufficient in order to handle discrete 
entropy inequalities. The main interest in such an estimate stems from the simplicity of its deriva- 
tion. 
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5.2 Discrete entropy estimates 

We first focus on the derivation of the discrete entropy inequalities and then the required weak 
estimate. The passage to the limit in the discrete inequalities is the subject of the following section. 
After Crandall and Majda [2T], assumptions (3.11 )-( 3.13) on the numerical flux functions g e ,K 
are known to yield a full set of discrete entropy inequalities for scalar conservation laws. Here and 
in the light of Section [2j the scalar conservation laws of concern have to be found locally at each 
edge e in 77i, and take the generic form 



d t w + V ■ f(w, v) = 0, 



(5.1) 



for a given v € 1R. Associated entropy pairs were defined earlier inJ2.40-(2.44). The inequalities 



stated below are naturally built from the subcell states w^~ e ' (4.1) of Lemma 
regard may be understood as subcell entropy inequalities 



4.1 



and in this 



Lemma 5.1 (Entropy inequalities per cell). Let (U, J-) : 1R —> IRx lR d be any convex entropy pair 
for the scalar conservation law (5.1), where e denotes any edge in dK for an arbitrarily K inTh- 
Then there exists a numerical entropy flux function G &i k • JR 2 ~> JR that satisfies the consistency 
property 

G etK (w,w;v e ) = F{w,v e ) ■ u K , e , (5.2) 



the conservation property 

G e ,K(w,w e ;v e ) = -G eiKe (w e ,w;v e ) 
for all reals w and w e , so that the following discrete entropy inequality holds 

U{w n £n-U{™ n K,e) + —^(Ge,K(w^, e ,W^, e ;V e )-T(w n Kie ,V e ) ■ V K , e ) 

UK,e KM V ' 



< 0. 



(5.3) 



(5.4) 



We refer the reader to [3T] for a proof of this classical result. As already claimed, the weak 
estimate will not allow to pass weakly to the limit in arbitrary numerical entropy flux-functions. 
We thus propose to merge inequalities (5.4 1 in such a way that solely exact entropy flux-functions 
J-(w, v e ) ■ VK,e enter the weak form. 

Lemma 5.2. Let (f> be any non-negative test function in T>{IR* + x IR d ). Define for any edge e in 
Th, the average 

1 



(x, i)dxdt. 



(5.5) 



Then, the following discrete weak inequality holds 

£ a ^( w «!- 1 '")- u{w n Kie ))<t> n e \K\- 



T E £-7 7 «,e,«e)-I'X, 8 #|e|<0. 
KeTh eedK 



(5.6) 



The proof is postponed to the end of this section. We shall easily deduce from the discrete 
inequality (5.6) the following continuous (in space) inequality. 



Proposition 5.3. The finite volume approximation (3.1J h )-(3.11) obeys at each time level t n the 
following (discrete in time) entropy inequality 



g £ a K ^U{wl + ^-)-U{w K ^y:\K\ 

Q{ul,v{x)) ■ V<p(x,t) + (f>(x,t)d v Q(u%,v(x)) : Vv(x)dxdt 



KeTu eedK 



(5.7) 



/]t™,t"+ 1 [xlR d 

< 0(h)r \\4>\\ w i,^ {]tn ^+i [xRd] \supp{4>)\. 
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Th e proof, at the end of this section, essentially makes use of the uniform sup-norm estimate 
(4.13) for the sequence together with the regularity assumption w n € W 2, °° . 

Clearly, the Young measure \x can tackle the weak limit of the space derivatives involved in 
inequality (5.7) extended to any time interval (0,T), T > 0. Such a claim then naturally rises the 
question of passing to the weak limit in the discrete time derivative. The latter is conveniently 
decomposed as 



£ Yl a K , e (u{w n £n-U{wl <e ))<t> n AK\ 



KeTh eedK 



' T h 



i<eT h 



KeTh eedK 



where 



W- £ Ya K , e (u(w n K +1 )-U(w n £>-))€\K\ 

KeT h eedK 

-U(w%))fi\K\, 



b K = E aK > e< < 
eedK 



(5.8) 



(5.9) 



The last two error terms entering the righ-hand side of (5.8) are devoted to sum up 



E E E oc K , e {u{ W ^)-U{wl^-))€\K\ 

n>0 KeT h eedK 



(5.10) 



n>0 KeTh eedK 



with other error terms in the right-hand side of the discrete entropy inequalities (5.7). The former 
must therefore be proved to go to zero with h. 



Lemma 5.4. For any polyhedron K in Th, one has 

a K , e (u{w^ e )-U(w n K ))<i>: < 0(ti< 



W 1 -°°Qt n ,t n + 1 [xK)i 



(5.11) 



eedK 



while 



eedK 

^- ( I n+1 n+l,-|2 

< -O-ul }^ a K,e\ w K ~ w 



J K, 



eedK 



(5.12) 



ra+1, 
' J K,e 



eedK 



|i=o(] t n it n+l[ x Jf), 



where o~u denotes some convexity-like modulus oflA: U"(u) > o~u > 0, for all u 6 (m, M) where 
the bounds m,M were introduced in (3.18) in agreement with the maximum principle \^.13 ). 



Proof of Lemma 5.2 Let e be any edge in Th and K, K e the associated pair of adjacent polyhe- 
dra. Multiplying the subcell entropy inequality (5.4) valid for K by a.K,e\K\ and the companion 
inequality for K e by atK,e\Ke\, we get 



OLK, 



m(u{w n £n-U{wX e j) + a Ke ,e\K e \(u(w n K + ^-)-U(w n Ke>e ) 

- T(T(w^ e ,V e ) ■ V K ,e + T{w n Kce ,V e ) ■ VK c ,e^j |e| < 



thanks to the conservation property (5.3) satisfied by the numerical entropy flux-functions. Mul- 
tiplying the above inequality by the discrete test function </>™ (5.5), then summing over the edges 
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J K B ,e) jfel-Kel 



e in dK and the polyhedra K in 7~h yields 

E E ^e^C^e 1 '") " W«, e ))#l*l 
+ E E a ^( W «te")-^ 

-T E E (-^(^Sr.ei u e) ' KK> + ffw^e^e) • ^K e ,ej0el e l < °- 
KeTh eedK 

To conclude the proof, we notice the following two identities 

E E ^e(W«+ 1 '-) -U{w^ e ))€\K\ 
KeT h eedK 

= E E a ^^( W «t,e")- W (^e,e))^el^e|, 

ifeTh ee9K 

and 

KeTh eedK KeTh eedK 



□ 



Proof of Proposition\5.3\ We begin with the discrete inequality (5.6 1 of Lemma 5.2 and specif- 



ically considerb the flux balance E E J~( w K,e' v e) ' v K,e<l>e\ e \- O ur purpose is to shift the 

KeTh eedK 

mathematical expressions under consideration from the w to the u variable. Hence let us write 
J(aj^ e ,D e ) = J : {w{u 7 j ( , v e ),v e ) = Q(u K ,v e ) with Q(u,v) the exact entropy flux introduced in 
(2.40), which we repeat component- wise as Qi(u, v) — f u U'(Co(9, v))dgCi(9, v)d0, 1 < i < d. We 
then recast the flux balance as 

E T ( W K,e^ e )-V K A n M 
dK 

Q«,v K )- £>e|eKe+ E (Q(u n K ,v e ) - Q(l 



eedK 



(5.13) 



U K, V K) ■ V K ,e0 e |e|, 



eedK 



eedK 



where the average of the states v e is defined by vk = &K,e v e- I n view of a representation 



eedK 



formula for V</> (similar to the one in Remark 3.1 derived for Vw), the average form (5.5) for 
yields 



E ^eWK^e = ~ ( E / <t>{x,t)vK,edx)dt = - / 



V(f>(x,t)dxdt, (5.14) 



eedK 



eedK 



so that, from (5.13), 



E ^Ke^e) ' VK,e4>7\ e \ 



eedK 



~ I [ Q( U K> V K ) ■ V4(X, t)dxdt + E (QKr, «e) ~ e(t&, «JC)1 ' ^,ede 



(5.15) 



The treatment of the last remaining discrete term relies on the following identity: 
Q(u n K ,v e ) - Q(u n K ,v K ) = I d v Q(u n K ,v K + s(v e - v K ))ds (v e - v K ) 
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which leads us to rewrite (5.151: 



i r tn+1 r 

V T{w n K e , v e ) ■ v K A n M - - / / Q{u K , v K ) ■ v<Kx, t)dxdt 

eedK T Jt " jK 

= d v Q(u K ,v K ) : <t> n e {ve-v K )®v Kte \ 



eedK 



„1 

+ {9vQ(u K ,v K +s(v e -v K ))-d v Q(u K ,v K ))ds^ : ((u e - v K ) <g) u K , e ) |e|. 

eedK 

(5.16) 

The matrix (v e — VK)®^K,e\^\ with size Lx d appears as a discrete representation for the continuous 
function Vu. The first term in the above right-hand side is rewritten as: 



d v Q(u K ,v K ) ■■ ( ^2 0e( Ve ~ VK "> ® v i<A e \) 

eedK 

= 4>k d v Q{u n K ,v K ) : ( (v e - v K ) ® ^if,e|e| 

+ (E (€-<f> K )dvQ(u K ,v K ) : ((« e -^)®^,e)|e|), 



(5.17) 



eedK 



where the discrete flux function <$>\ is obtained by averaging: ip K = a^^". On one hand, 

eedK 

owing to the identity (u e — vk) ® VK,e\&\ — u e ® v K,e\ e \ we get 



eedK 



eedK 



( 1 

floCCuSr.wjc) : (5^ (« e - ® fjc,e|e|) = d v Q(u K , v K ) : ( - J Vv(x)dtd 



(5.18) 



again thanks to the representation formula in Remark 3.1 (for Vi>). On the other hand, the latter 
error term in ( |5.17[ ) is described by 



| (€-r K )d v Q(u n K ,v K ) : ((v e - v K ) ® i^,e)|e| 

<0(1) SUp |(^-05-)(«e-«ic)|)pjC 

< 0(/&)||V0|| i ~ at », t »+i [xJ r ) p* < O(fcir) ||V0|U» (]t „ t „ +1[xK) |X|. 



(5.19) 



Here, we have successively used the sup-norm estimate (4.13) satisfied by Uh, the definition of the 
perimeter px of K, the estimate 



\v e -V K \< E OL K ,e\ v e ~ V e >\ < 0(h K ) 



(5.20) 



e'edK 



from the definition of Vk and the regularity property vq £ W 2,co , a similar estimate |0™ — </>JJ| < 
OQik) and finally the non degeneracy assumption (3.1 1 on the triangulation 7/t- Involving (5.181- 
(5.19), the identity (5.17) yields the following estimate 



1 r* r 

d v Q(u K ,v K ) : 4>2{v B -v K )»v KiB \e\ \ -- / / ^-^Q«,^):V« 



(x)dtdx 



<0(h K ) ||V</»|| L oo (r)t „ +1[xX) |^|. 



(5.21) 
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For the final error term in the flux balance (5.16), we have the following bounds: 

€ I 

eedK 



^2^e {dvQ(uK,v K + s(v e -v K ))-d v Q(u%,VK))ds:((y e -v K )®v K<e )\e\ 
edK Jo 

< 0(1) Sup \v e - Vif| 2 (pK||0||L~(]t".t"+i[xK) J 



eedK 

<0{h K )\\<l>\\L°<>{]t*,t^lxK)\K 



where we have used the regularity of the entropy flux Q, the sup-norm estimate (4.13 ), the estimate 
(5.20) satisfied by \v e — vk\, and the non-degeneracy assumption (3.1) on the triangulation Th- 
To summarize, we have obtained the estimate for the flux balance on a single cell: 



t" jk 



(Q{u n K ,v K ) ■ + <t>KdvQ{u 1 k,v K ) : Vv{x)dtdx 



eedK 

<0(h)U\\ W ^ (]t n, t n +1[xK) \K\. 



(5.22) 



From the discrete weak entropy inequality (5.6) we recall that 



Y / a K ,e(u(w n K + e 1 >-)~U(w^ e ))^\K\-T ]T J2nw n K!e ,v e )-p K>e ^\e\<0, 



eedK KeT h eedK 

the sum of ( |5.22[ ) over all cells K on the triangulation Th gives 

>K,e))€\K\ 



J2 a ^( u( - w Z 1 n-^^ 

KeTh eedK 



r r 

/ ( J2 Q{u n K ,v K )-V<t> + <j>d v Q{u K ,v K )-.Vvdx) 
Jtn V e r„ Jk 



dt 



<0(h)r Y \\<l>\\wi.°°Qt»,t"+i[xK)\K\ < C( /l ) T ll0llwi.~(]t",t"+i[xiR d )l su PP( < / , )l- 
Ken 



□ 



Proof of Lemma 5.4 We first establish the estimate (5.11 1 and consider the following decomposi- 
tion involving again the {aR,e}{e,eedK }-average <p K of the 0™ (5.9): 

X>*,e("«, )-W(«&))# 
eedK 

= Y. a KJu{w n K ^)-U{w n K ) 



eedK 



n \ \ ( in in \ 

" , -9k) 



n K ( K Y. a ^ U ^K,e)-U{w K ) 



eedK 



from which we deduce the following bound: 



eedK 



<0(h K )\\V<i>\\L°°Qt»,tn+ilxK) SUp \W% -W%\ 



eedK 



(5.23) 



O 



(1) \\9\\L~( ]t n, t ^ [x K)( a K , e U(w^ e )-U(w n K ))., 



eedK 
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in view of the sup-norm estimate (4.131 satisfied by Uh, the estimate |</>" — <fi K \ < 0{Kk) and the 
convexity of the entropy U(w). The hrst error term in (5.23) is given the following bound: 



\Wk iS ~ W%\ < ^ aK,e'\C (uK,Ve')-Co{u%,V e )\ 
e'edK 

< 0(1) sup |tv - v e \ < 0(h K ), 

e'edK 



while the second one may be handled as follows: 

]T a K , e U(w n Kte ) -U(wl) =U'(w n K )( £ a K> , 



(5.24) 



eedK 



e™K,e ~ W K 



eedK 



+ J2 fu"{wl, e + S (w K - W Kte ))d S (w Kte - W n K f 

JO 



(5.25) 



eedK 

<0(1) sup K, e -<| 2 < 0(h* K ), 

eedK 



in view of (3.15) w K = J2 e edK aK - eW K e an< ^ ^he estimate (5.24). Gathering bounds (5.24) and 



(5.25) yield the expected estimate (5.11) in Lemma 5.4 



We now derive the companion estimate (5.12), by starting from the decomposition 

^ aK , e (w« +1 )-w« v 



eedK 



b n K (u{w n K +i ) £ aj r , B w«+ 1,_ )) + E «x, e (w« +i ) -u(w n £n) (€ - r K ), 



eedK 

and observing, on one hand, 



eedK 



£ ttJf , e (wK +i )-w«'i) (c - r K 



eedK 

<0(l) ]T «K,e|C-^IKTe , "-^" +l1 



A" 



< 



^(^(E^Kle 1 ' 



n+l 
K 



eedK 



|L°°(]t",t"+ 1 [x A") 



and, on the other hand, 

eedK 

eedK 

+ E a ^ e [^"(w^- +s(w 



eedK 



Finally, in view of the convex decomposition (4.2) stating = J2eedK a K, 



n+l, 
z W K,e 



U( W K +1 ) ~ E a K,eU{w K + ,e' ) ^ ~ a U E a KA W K + ,e' 



W 



n+l|2 
K I i 



eS9A 



eedK 



where ou denotes the convexity like-modulus of U introduced in Lemma |5.4| This concludes the 
proof. □ 
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5.3 Entropy dissipation rate and strong convergence 



The proposed estimates obtained in Lemma |5.4| deserve a few comments. Plugging first estimate 
( 5.1 f [ ) in ( |5.10[ ) will be easily seen to yield the following upper-bound 



E E E 

n>0 KeTh eedK 



s (w« ie ) -U(w n K ))^\K\ < O(/i)||0|| wl ,oo (K+xJRd) |supp(0)| 



that obviously suffices to conclude. By contrast and turning considering (5.121, a crude upper- 
bound based on the sup-norm estimate (4.14), say 



E aK << 

eedK 



u{w n K +l )-u{w n K + ^-))€. < o(hm\ w ^ Qtn<t , 



would result in the useless estimate 

E E E u Kte (u{wl +1 )-U{wl^n)&\K\ <O(l)||0|| w i.» (B+xB , ) |supp(0)|. 

n>0 KeTh eedK 



Proving that the err or ter m of concern in (5.10) actually vanishes with h requires therefore in turn 
a sharper control in (5.12 ) of the oscillations of the around their mean value u^ +1 . Such a 

control over these discrete oscillations results from a sharp evaluation of the discrete entropy rate 
of dissipation. 

Proposition 5.5. Let T > be any fixed time and let Nt £ N be the floor of T/t we denote 
[T/t]. Then, for any (time independent) non negative test function ip € T>(lR d ), the finite volume 
approximation (3.14)-(3.17) obeys the following estimate on the discrete oscillations: 



E a ^K +1 



n=0 KeT h eedK 



W 



n+1,- |2 
Ka 



iPk\K\<0(1), 



(5.26) 



where ipx reads ipK = /J aK,eipe 

eedK 



i>e 



ip(x)da 



Equipped with (5.26) we obtain the following entropy dissipation rate. 
Corollary 5.6. The sequence Uh satisfy the entropy like inequality 

U(C (u h , v))d t (j>(x, t) + Q(u h , v) ■ + <j)d v Q{u h , v) : Wvdxdt > 0(h^ 2 ), 

M+xR d 



(5.27) 



for any (smooth) convex entropy pair (U, Q) : 1R — > 1R x lR d introduced in ( |2.42 1 and (2.40). 

Equipped with the above inequality valid for any entropy pair (U,Q), we easily deduce that 
the Young measure /i = fi t ,x associated with the sequence (v,h)h>o is an entropy satisfying measure 
valued solution. In other words the uniformly bounded L°° sequence (uh)h>0: as announced at 
the beginning of this section, it is easy to check that the inequation ( |5.27[ ) becomes as h tends to 
the following inequation satisfied in the weak sense: 



d t {tiM{C Q (;v))) + V^/i, fi(.,tO) - (fi,d v Q(-,v)):Vv < 0. 



(5.28) 



Relying on a direct extension of DiPerna's uniqueness theorem [15] . we can deduce that the 
entropy measure-valued solution fj, t , x reduces to a Dirac measure 5 u ( t , x ) concentrated on a function 
u = u(t, x) since the initial data /io coincides with the Dirac measure S Uo (where uq is the initial 
data in the Cauchy problem (2.33)). Proving that the inital data uq is correctly handled amounts 
to show that for every compact subset K, of M we have 



lira 

t->0+ 



(fJ, s , x ,\id — u (x)\) dxds = 0. 



(5.29) 



o Jic 



28 



The condition ( 5.29 1-( 5.28 1 reduces to a Dirac measure concentrated at u(t,x), the Kruzkov en- 
tropy solution of ( 2.33 )-( 2.42 1 with same initial data u$. In other words, for all time T > 
and for all compact JC in M, the scheme converges strongly in Lf ((0,T) x JC) to the solution u. 
Theorem 3.4 of this paper is thus now established. 



Proof of Proposition \5.5\ We start from the discrete in time weak formulation (5.7 1 stated in 
Proposition |5.3| 

Y Y a K , e (u(w n K + e x n-U(w%,e))€\K\ 

SK,»W) • V<f>(x,t) + 4>(x,t)d v Q(u%,v(x)):Vv(x)dxdt 



KeTh eedK 



/]t",t"+ 1 [xiR d 

< 0(h)r II^H^i.ooQtn^n+i^^lsupp^)!, 



in which we plug the decomposition (5.8)-(5.9). A discrete test function ipx given for any given 
time-independent test function ip G T>(M d ) is considered. We then get 

Y (u{w n K +1 )-U{wl))i, K \K\ 

QK, v[x)) ■ Vip(x) + ip(x)d v Q(ul,v(x)):Vv(x)dxdt 



< y Y. a '<^ w K+ 1 )-u{w^-))^\K\+ y Y. a ^(u{w^ e )-u{w K ))y e \K\ 

KeT h eedK KeTh eedK 

+ OWrH^H^i.oo^lsupp^)!. 
Invoquing estimates (5.11 )- |5~l~2 | then yields 

Y {u{wl+ l )-U{w^-)y K \K\+a u Y J2 a ^\ w K +1 - w Z X '~\ 2 ^K\K\ 

KeTh KeTh eedK 

<0(fc)r||^| W r ll » (Rd) |Bupp(V)|+ 0(h) ]T \\VtP\\l~(k)\K\+ 0(h 2 ) Y \\nw^(K)\K\ 

KeT h Ken 



Q(v%,v(x)) ■ Vip(x)+i/>(x)d v Q(u%,v(x)):Vv{x)dxdt. 



Observe that due to the estimate (4.13), the last contribution in the above right-hand side can be 

Uk\K\ 



given the following crude estimate 0(r 



jyi.oo^d). Henceforth, we deduce that 



e\w K 



n+l ,,,n+l,-|2 



10 



K.i 



KeTh 

< o(h)U\\ 



KeTh eedK 



Summing over time indices n £E [0,Nx] with iVy = [T/t] for a fixed time T > 0, we get 



f NT ( 

i U(w h (x,T))tp h (x)dx + auY ( z2 aic 

IR d ^ — n TSr-T -rats- 



e\W K 



n+l ,„n+l,-|2 



W 



K.i 



< 



n=0 KeTh eedK 

U(w (x))ip h (x)dx + 0{l)T\\il}\\ w i,^ {]Rd) , 



which is the required result. 



□ 
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Proof of Corollary \5.6\ We start from ( 5.7 )-( 5.8 )-< 5.9 ) and consider the following discrete in time 
weak formulation for the time dependent test function (f> £ V(lRf x IR d ) and its discrete repre- 
sentation <f>\ 



(u(w- +1 )-U(w n K ))^\K\ 



K<£T h 

Q{ul, v(x)) ■ Vcj>{x, t) + (f>(x, t)d v Q{ul, v(x)) : Vv(x)dxdt 

l]t",t™+ 1 [x]R d 

< C , (/i)r||0|| w i, M(]t „ it „+i [xiRt!) |supp(</-)| + 0{h 2 ) \\<f>\\wi.°°Qtn,t»+i{xK)\K\ 

Ken 

+ o{h) E a ^i< +1 -<e 1,_ i m\\^ Qtn , tn+HxK) \K\. 

KeT h eedK 



where we have used estimates (5.11)-(5.12). Summing this inequality over time indices gives 



E E 

n>0 K€T h 



h n+1 — (h n 



Q{ul, v(x)) ■ V4>(x, t) + <P(x, t)d v Q{ul, v(x)):Vv(x)dxdt 

l]R+xR d [0.6V) 

<0(h)\\<t>\\ 

+ o a)E E (E a ^i^ +1 -^"i^n v ^i^(]*"^ + MxK)i^) 5 

n>0 KeTh eedK 

making use of the characteristic function \4> 01 Uo<t<T SU PP0K"5 a compact subset of M d , 
where T is a finite time such that supp(0(-,t)) = for t>T. Cauchy-Schwarz's inequality then 
yields the following crude upper bound for the last term: 



E E (. E <>- '<■;;• ; - 1 



w K^e' I X*)\\V<l>\\L°°(]t«,t«+i[x:K)\K\'r 

<(EE (E^i^ +1 "-^ix.) 2 Wr) 1/2 (E E 11^ 

n>0K£T h eedK 

<o(i)(E E (E^K 1 -*"! 1 )*; 



i>0 KeT h eedK 



li~(]t",t™+ 1 [xX) 



\K\t 



n>0 KeT h eedK 



n>0KeTh 
1/2 



1/2 



as a consequence of the convexity property of the otK,e~ average. The estimate (5.26) then yields 
with ip = X4> 

E E (E a ^>^ 1 - w K , "U«)l|V0|U- a tn^ + i [xX) |ir|r< 

n>0 fi-eTk eS9if 



Then routine arguments give the conclusion from (5.30|. 



□ 



6 Numerical experiments 

6.1 A two domain coupling problem 

In this first test, we consider an heterogeneous medium which occupies the spatial domain [—1, l] 2 
and is constituted by an annular inclusion T>i centered at the origin (0, 0) with external radius \/0.2 
and with internal radius VO.l, and by its complement set T>q. In these two domains, the following 
respective flux-functions are considered in term of the scalar unknown w = w(t,x): 

ff ,™ 2 [l\ ft , (w-0.9)* fl\ 
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The regularized color function v plotted in Figure [3b| provides us with a regularized version of the 
characteristic function of the domain T>\. The coupling condition between T> and T>\ takes here 
the form 

2w-(t, x) = w + (t, x), x G &D±, 

where w±(t, x) = limg_>o+ w(t, x ± 0v x ) and v x the exterior unit normal at x G &D±. 
The initial data plotted in Figure [3a| is piecewise constant: 



wo(x,y) 



1, x < -0.8, 
0, x > -0.8. 



The computations are performed on a Cartesian grid with 100 x 100 meshes, and the CFL number 
is chosen to be 0.5. 

In an homogeneous domain with the sole flux /o, such an initial data would develop a shock front 
moving with the speed vector 0.5(1, 1) T . In the present heterogeneous domain, this shock front 



has the same behavior only until it reaches the interface between both domains (see Figures 4a). 
The coupling condition at this interface is such that the value w = 2 arises then inside the domain 
T>+. In this second domain, where the flux under consideration is /i, we observe then a (curved) 
shock wave connecting the states w = 2 and w = and moving at the fixed speed given by the 
Rankine-Hugoniot relation, that is, 0.605(1, 1) T (see Figures 
goes outside the whole domain [— 1, l] 2 (see Figure 4g 



4c 



In Figures Wb 



a nd |4e|) . Finally, the shock front 
4fj and 4h | , we plot the 



4d 



it-variable, which is found to remain constant at each interface, as expected by the theory. 





(a) Initial data wq (b) Color function v 

Figure 3: Initial data for the multidimensional test. 



6.2 A three domain coupling problem 

In this second test, we consider three different domains, as represented by the two components of 
v (see Figure [5]). The domain T>2 is a triangular inclusion and the domain T>i is the complement 
of T>2 relative to an annular inclusion. The flux-functions under consideration are now 



/oH 



w 
~2 



AH 



w 
~2 



0.5 




/2H 



w 2 f0 



(6.1) 



and the coupling relations are given by the change of unknown (2.23) with 



9 (w) = w, 9 1 {w)=w/2, 9 2 (w)=w/3. 



(6.2) 
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We consider the same initial data as previously and, thus, we expect the state w — 2 to appear in 
T>\ and the state w = 3 in T> 2 . The results are represented in Figures [6a] to [6f] for successive time 
steps. Once again, the limiting solution as the time grows satisfies the expected coupling relation. 
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(a) Solution w at t = 0.5 



(b) Solution u at t = 0.5 





(g) Solution w at it = 4.5 (h) Solution u at t = 4.5 



Figure 4: Evolution of the solution fof^ifferent times : w (left) and u (right). 




1-8. 8-8. 6-8. +8. 2 0.2 8.4 0,6 8.8 1 1 -8. 8-0. 6-0. +8. 2 8 8.28.48.68.8 1 

(a) Domain T>\. (b) Domain T> 2 . 

Figure 5: Geometry of the three domains. 
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(a) Solution w at t = 1.0 



(b ) Solution w at t = 2.0 




(e) Solution watt — 5.0 (f) Solution w at t = 6.0 

Figure 6: Three domain evolution. Solution w. 
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